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1.  INTRODUCTION 


With  the  widespread  use  of  numerical  (hydro-)  codes  to  describe  the  behavior  of  bodies  impacting 
at  high  velocity,  it  is  vital  that  code  developers  and  users  be  familiar  with  the  various  numerical 
models  in  the  codes.  All  too  often,  effort  is  concentrated  on  the  "numerics"  of  a  code,  such  as 
integration  and  discretization  techniques,  artificial  viscosity  formulations,  mass  lumping  and  pressure 
averaging  techniques,  and  mass  advection  schemes,  with  the  same  rigor  being  absent  when  considering 
the  "physics”  of  the  code. 

Because  code  development  involves  significant  effort,  it  is  common  for  codes  to  evolve  or  be 
passed  down  directly,  rather  than  be  written  from  the  ground  up.  With  this  sort  of  evolution,  it  is 
possible  to  avoid  "reinventing  the  wheel."  Unfortunately,  this  mentality,  if  taken  to  the  extreme, 
precludes  the  ability  to  redesign  the  “wheel"  and  can  result  in  an  incompatible  combination  of  "wheel" 
and  code  "vehicle."  Common  problems  associated  with  material  model  evolution  are  enumerated 
below: 

(1)  Many  materia'  models  arc  based  on  assumptions  and  limitations  which  are  often  not  passed  on 
w  ith  the  model.  Thus,  material  models  which  arc  valid  only  for  a  particular  class  of  problems 
ic.g.,  small  compression,  ductile  behavior,  etc.)  can  be  erroneously  applied  to  more  general 
classes  of  problems. 

(2)  Goodness  of  a  material  model  is  often  judged  on  its  ability  to  solve  a  particular  set  of 
problems,  rather  than  on  adhcrancc  to  fundamental  principles.  Such  heuristic  approaches 
(e  g.,  use  of  non-associatcd  flow  rules)  can  be  exercised  with  caution  but  often  become 
institutionalized  and  taken  as  fact  by  an  unknowing  community. 

13)  Many  new  material  models  arc  derived  by  applying  ad  hoc  modifications  to  older  models, 
rather  than  rcdcriving  the  models  from  basic  principles.  When  this  is  done  carelessly, 
conflicting  assumptions  or  limitations  can  render  such  models  invalid. 

(4j  Material  libraries  tend  to  take  on  a  life  of  their  own  and  arc  passed  down  from  generation  to 
generation  without  the  accompanying  experimental  data  used  to  generate  those  libraries. 
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When  this  happens,  material  data  are  often  used  under  conditions  (e.g.,  strain  rates,  pressures, 
etc.)  where  those  data  arc  not  valid. 

One  aspect  of  material  modeling  which  is  extremely  important  to  impact  codes,  and  to  which  all  of 
the  abovementioned  problems  arc  relevant,  is  the  equation  of  state  (EOS).  In  this  report,  equations  of 
state  for  use  in  impact  codes  are  discussed  in  general  terms,  with  the  Mie-Grilneisen  equation  of  state 
used  for  specific  illustrative  purposes  since  it  is  the  workhorse  EOS  of  the  current  generation  of 
hydrodynamic  wave  codes.  Basic  derivations,  assumptions,  and  limitations  of  applicability  arc  all 
discussed,  with  the  prime  focus  being  on  cquation-of-statc  stability. 


2.  THERMODYNAMIC  STATES  AND  THE  EQUATION  OF  STATE 

To  completely  describe  the  volumetric  behavior  of  a  material,  the  values  of  the  thermodynamic 
state  variables  (i.c.,  the  thermodynamic  coordinates)  arc  needed  for  every  state  of  the  material. 
Thermodynamic  state  variables  arc  properties  which  arc  only  dependent  on  the  state  of  a  material  and 
not  on  the  path  taken  to  arrive  at  that  state.  Examples  of  typical  thermodynamic  state  variables 
include  pressure  and  temperature,  as  well  as  specific  volume,  internal  energy,  enthalpy,  entropy,  etc. 

To  describe  the  state  of  a  simple  substance,  excluding  such  phenomena  as  phase  change, 
dissociation,  or  electromagnetic  work,  a  knowledge  of  the  material’s  pressure,  as  a  function  of  density 
and  temperature,  is  sufficient  to  completely  describe  the  state  of  a  material.  In  general,  these 
thermodynamic  state  data  are  related  in  one  of  three  forms. 

Tabulated  experimental  data  is  one  popular  form  of  representing  thermodynamic  state  data.  Steam 
tables  arc  probably  the  most  widely  used  form  of  tabulated  thermodynamic  data.  Data  are  tabulated 
over  a  large  range  of  states,  and  interpolation  becomes  necessary  to  obtain  data  between  the  tabulated 
state  points. 

Of  a  similar  nature  to  tabulated  thermodynamic  data  is  graphical  thermodynamic  data.  The 
Mother  diagram,  which  graphically  represents  the  states  of  water  in  its  various  phases,  is  the  classic 
example  of  graphically  represented  thermodynamic  data.  To  the  experienced  user,  graphical  state 
tables  provide  a  fast,  accurate  means  of  evaluating  the  change  in  state  variables,  when  something  is 


2 


known  about  the  thermodynamic  process  in  question.  Like  tabulated  data,  extensive  experimental 
testing  is  necessary  to  generate  »be  data  found  in  graphical  thermodynamic  state  tables. 

The  third,  and  possibly  most  common,  way  to  represent  thermodynamic  state  data  is  in  analytical 
form,  and  is  referred  to  as  an  equation  of  state.  Traditionally,  the  EOS  expresses  pressure  in  terms  of 
density  and  temperature.  Describing  material  by  way  of  an  EOS  has  many  advantages.  The  form  is 
compact,  much  more  so  than  exhaustive  tabulated  data.  An  EOS  may  be  mathematically  codified,  thus 
making  it  the  most  attractive  on  the  basis  of  programming  considerations.  Additionally,  the 
interrelationship  of  state  variables  may  be  mathematically  derived  from  the  EOS,  thus  providing  a 
more  fundamental  understanding  of  the  thermodynamic  phenomena  at  work. 

On  the  other  hand,  equations  of  slate  can  suffer  from  severe  deficiencies,  which  can  make  their  use 
prone  to  error.  Invariably,  equations  of  state  arc  formulated,  with  a  host  of  inherent  assumptions  about 
the  nature  of  the  material  behavior.  Common  examples  of  such  assumptions  include  constant  specific 
heats,  ideal  gas,  etc.  While  these  assumptions  may  be  valid  over  a  small  domain  of  possible  states, 
they  are  often  poor  at  describing  the  material  under  extreme  thermodynamic  conditions  like  high/low 
pressure,  phase  change,  and  the  like.  This  limitation,  in  and  of  itself,  is  not  bad.  However,  those  who 
use  these  equations  arc  not  always  cognizant  of  the  conditions  under  which  a  given  EOS  is  valid. 

Thus,  equations  of  state  may  be  used  in  a  thermodynamic  domain  where  the  EOS  parameters  are  no 
longer  valid.  Also,  there  can  be  a  problem  with  EOS  forms  which,  though  capturing  the  essence  of 
the  material  behavior,  suffer  from  "higher  order"  inconsistencies  that  only  manifest  themselves 
obliquely. 

When  it  comes  to  hydrocode  modeling  of  material  deformation,  where  the  EOS  is  just  a  means  to 
a  computational  end,  sophistocatcd  equations  of  state,  though  available  in  the  literature,  arc  rarely  used 
because  of  simplicity  (i.c.,  efficiency)  considerations.  This  focus  on  computational  simplicity  merely 
serves  to  enhance  the  likelihood  of  EOS  failure  in  hydrocode  compulations. 

All  equations  of  state  need  data  with  which  to  drive  themselves.  In  the  case  of  the  Mic-Griinciscn 
EOS,  as  used  in  today's  hydrocodes,  thermodynamic  state  variables  arc  expressed  relative  to  those 
states  found  along  an  experimentally  determined  pressure-volume  curve  which,  for  the  case  of  impact 
data,  is  usually  chosen  as  the  Kugoniot  curve.  Thus,  before  the  Mie-Griineisen  EOS  can  be  discussca 
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in  detail,  the  reader  must  acquire  an  appreciation  of  the  Hugoniot  and  the  thermodynamic  process  of 
shock  transition. 


3.  SHOCK  TRANSITION:  A  THERMODYNAMIC  PROCESS 

A  shock  is  a  volumetric  disturbance  which  travels  faster,  in  a  medium,  than  the  bulk  speed  of 
sound  (i.c.,  sonic  velocity)  characteristic  to  the  state  of  that  medium.  A  shock  is  characterized  by  a 
shock  front,  which  itself  is  very  thin  (thickness  on  the  order  of  micrometers),  and  is  analytically 
idealized  as  a  discontinuity  in  physical  properties  such  as  density,  pressure,  and  energy.  The  thinness 
of  the  shock  coupled  with  the  high  velocity  of  the  shock  front  ensure  the  adiabatic  nature  of  shock 
transition.  Large  amplitude  disturbances  arc  not  necessarily  shocks  if  the  thickness  of  the  disturbance 
is  relatively  large.  Correspondingly,  the  amplitude  of  a  shock  need  not  be  necessarily  large  as  long  as 
the  disturbance  is  traveling  faster  than  the  sonic  velocity. 

Shock  transition  may  be  thought  of  as  an  adiabatic  irreversible  thermodynamic  process  which 
permits  the  state  of  a  material  to  be  altered  in  a  fashion  similar  to  iscntropic,  isothermal,  isochoric, 
iscnthalpic,  and  isobaric  processes.  Whereas  all  of  these  other  thermodynamic  processes  arc 
characterized  by  the  constancy  of  a  particular  state  variable,  the  process  of  shock  transition  is 
characterized  only  by  the  existence  of  a  shock  wave.  No  thermodynamic  state  variables  remain 
constant  during  shock  transition.  Nonetheless,  the  shock  transition  process  can  be  as  well 
characterized  as  any  of  the  other  thermodynamic  processes,  as  will  be  shown  from  the  following 
derivation  of  the  basic  shock  relations. 


4.  BASIC  SHOCK  EQUATIONS 

Consider  a  shock  disturbance  traveling  at  velocity  U,  into  a  stationary  medium,  with  density, 
pressure,  and  specific  internal  energy  of  p0,  p0,  and  E0,  respectively.  Behind  the  shock,  the  material 
has  acquired  a  panicle  velocity,  given  by  up,  traveling  in  the  same  direction  as  the  shock  disturbance. 
The  density,  pressure,  and  specific  internal  energy  arc  p,  p,  and  E,  respectively.  This  situation  is 
depicted  in  the  Laboratory  Frame  diagram  of  Figure  1.  If  an  observer  could  be  situated  right  on  the 
shock  wave,  so  that  the  shock  were  to  appear  stationary,  the  situation  would  appear  like  that  depicted 
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Figure  1.  Deniction  of  i -ahoratorv  and  Shock  Reference  Frames,  Used  to  Derive  Continuity 
Momentum,  and  Energy  Relations  for  Shock  Transition  (Zuckcr  1977). 


in  the  Stationary  Shock  Frame  diagram.  The  continuity,  momentum,  and  energy  relations  are  derived 
below  by  considering  an  infinitesimally  thin  conlol  volume,  which  fluxes  mass  dM,  over  a  time 
increment  dt,  which  straddles  the  shock  front,  in  the  stationary  shock  frame.  The  control  volume 
projects  unit  area  A  onto  the  shock  front. 

CONTINUITY: 

Mass  In  :  pa  Us  A  dt 
Mass  Out:  p  (Ua  -  u ^  A  dt 


Or,  in  terms  of  compression  p: 


where  p  =  —  -  1 
P0 


MOMENTUM: 

Dccelcrativc  Impulse:  (p  -  pj  A  dt 

Momentum  Decrease:  dM  [Us  -  (1/,-u^y  =  fp,  U,  A  dt)  [Ut  -  (Ut-u^J 

P  ■  Pa  =  P0  V,  V  f°r  Po  =  P  =  P.  Us  Up 

ENERGY: 

Energy  In: 

Internal  Energy  In  :  Eg  dM  =  Eg  ( p0  U,  A  dt) 

Kinetic  Energy  In  :  1 12  dM  U}  =  112  (pt  Ul  A  dt)  Ut2 

Net  Work  on  Control  Volume:  I F  dx  =  (p0  A)  (U,  dt)  -  (p  A)  [(U^u^  dt] 

Energy  Out: 

Internal  Energy  Out  :  E  dM  a  E  (p#  U,  A  dt) 
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Kinetic  Energy  Out  :  1/2  dM  u2  =  1/2  (p0  U,  A  dt)  ( U-u J2 

Internal  Energy  Accumulated:  0 


E  -  Ea  =  1/2  (p  +  pJ(l/Po  ■  1/p);  for  p0  =  0,  E-  Ea  =  1/2  p(llp0  -  Up) 


Or,  in  terms  of  compression  | a: 


E  -  Eo 


(P  +Pe)  P  . 
2po(l  +  Ii)  ; 


forpo~0,  E-E0 


P  H 

2pa  (1  +  p) 


5.  THE  IIUGONIOT 

Tlic  momentum  relation  derived  above  expresses  pressure  in  terms  of  shock  and  particle  velocities, 
Us  and  up,  respectively.  If  the  nature  of  the  shock  or  particle  velocity  relationships  are  explicitly 
known,  then  the  pressure  can  be  expressed  in  terms  of  specific  volume  v,  density  p,  or  alternately, 
compression  p.  This  pressure-volume  relationship  for  the  shock  compression  process  is  called  the 
Rankinc-Hugoniot  equation — or  more  simply,  the  Hugoniot — and  is  the  locus  of  states  which  can  be 
achieved  through  the  shock  transition  of  a  material  from  a  given  (pQ,  ve)  state. 

It  is  vital  to  note  that  the  Hugoniot  docs  not  represent  the  actual  path  of  states  through  which  a 
material  progresses,  when  transitioning  from  the  (p0,  vj  state  to  (p,  v),  but  rather  the  locus  of  final 
(p,  v)  states  which  can  be  achieved  through  shock  transition.  Also,  since  we  know  shock  transition  to 
be  a  compressive  procccss  only,  a  particular  Hugoniot  curve  has  one  endpoint:  the  pre-shocked 
(p0,  vG)  state.  From  this  endpoint,  the  Hugoniot  proceeds  in  the  direction  of  increasing  pressure  and 
decreasing  volume. 

Examining  the  shock  energy  relationship  for  the  ease  of  an  infinitesimal  shock,  it  is  seen  that 
dE  =  -pc  dv,  which  exactly  equals  the  reversible  work  that  is  done  on  the  material  during  the 
infinitesimal  compression.  This  condition  makes  the  process  reversible  adiabatic,  thus  iscntropic.  The 
implication  is  that  at  its  endpoint,  the  Hugoniot  is  tangent  to  the  isentrope  through  that  point.  Except 


7 


at  its  endpoint  however,  the  Hugoniot  slope  (with  respect  to  density)  is  necessarily  sleeper  than  the 
isentrope  through  a  given  point. 

Unlike  other  thermodynamic  processes,  for  example  iscntropic  compression,  where  knowledge  of  a 
single  state  on  an  iscntropic  (p,  v)  curve  is  sufficient  to  completely  describe  that  iscntropic  curve, 
knowledge  of  a  state  point  on  a  Hugoniot  curve  is  not  sufficient  to  define  the  Hugoniot  curve.  That  is 
to  say,  there  arc  many  (p0,  vj  states  which  can  possibly  shock  transition  to  state  (p,  v).  This  is 
identical  to  saying  that  there  arc  many  Hugoniots  through  any  state  point,  as  we  shall  prove  in  the  next 
section. 


6.  SOME  PROPERTIES  OF  THE  HUGONIOT 

To  prove  the  assertion  that  there  is  no  unique  Hugoniot  through  the  state  point  (p,  v),  let  us 
assume  the  contrary  and  show  an  inconsistency.  Consider  the  Hugoniot  referenced  to  state  (pQ,  v0) 

(we  will  hereafter  call  this  the  (p0,  vj  Hugoniot)  which  goes  through  stale  point  (p,  v).  Consider  also 
another  Hugoniot,  referenced  to  state  (p,,  v,),  which  also  goes  through  state  (p,  v).  If  one  assumes 
that  there  is  a  unique  Hugoniot  which  goes  through  the  state  (p,  v),  then  one  is  forced  to  conclude  that 
both  states  (p0,  vQ)  and  (p,.  v,)  lie  along  this  unique  Hugoniot.  If  (p,f  v,)  is  in  fact  on  the  (p0,  vo) 
Hugoniot,  then  shock  transition  can  theoretically  occur  from  the  state  (pc,  vj  to  (pj,  v,). 

To  show  the  inconsistency  which  arises  from  this  assumption,  consider  the  changes  in  specific 
internal  energies  which  arise  when  shock  transitioning  from  states  (p0,  vc)  to  (p,  v),  from  (p,,  v,)  to 
(p,  v)  and  also  from  (p0,  v0)  to  (p,f  v,). 

E  -E0-  1 !2  (pa  +  p)  (v0  -  v)  , 

E  *  E,  =  1/2  (p,  +  p)  (v,  -  v)  , 

E,  ■  E„  =  / 12  (pQ  +  p,)  (v0  -  v,)  . 


Subtracting  the  second  two  cquatioas  from  the  first  should  produce  an  identity,  but  it  can  be  readily 
seen  that  not  all  of  the  terms  cancel:  * 
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0  =  1 12  (pa  (v,  -  v)  +  p  (va  -  V,)  -p,  (Va  -  y)}  . 

This  can  be  reduced  to: 

(Pi  -  Pc)  1  <vo-  vi)  ~  (P  ‘  Pi)  1  fvi  '  v )  ■ 


Not  only  is  this  expression  not  an  identity,  but  it  will  only  reduce  to  a  true  statement  if  the  Hugoniot 
coincides  with  a  Rayleigh  line,  which  is  line  of  constant  slope  in  the  (p,  v)  plane.  Additionally, 
thermodynamics  tells  us  that  there  exists  a  local  maxima  of  entropy  along  a  Rayleigh  line,  where  it  is 
tangent  to  an  iscntropc.  Thus,  the  Hugoniot  cannot  possibly  lie  along  a  Rayleigh  line,  for  to  do  so 
would  imply  decreasing  entropy  with  increasing  shock  strength  as  the  shocked  state  moved  beyond  the 
local  entropy  maxima  on  the  Rayleigh  line. 

Thus,  the  inconsistency  has  been  shown,  and  one  can  conclude  that  the  assumption  that  both 
(p0,  vc)  and  (p,,  v,)  lie  on  the  same  unique  Hugoniot  must  be  false,  even  though  the  Hugoniots 
originating  at  both  these  points  intersect  at  the  point  (p,  v). 

A  corollary  to  this  point  is  that,  if  one  were  to  proceed  from  state  (p0,  vj  towards  a  specific 
volume  of  v,  by  way  of  two  shocks,  the  final  pressure  will  not  be  the  same  as  when  specific 
volume  v  is  achieved  through  a  single  shock,  since  the  new  Hugoniot  originating  at  the  intermediate 
shocked  state  will  differ  from  original  (p0,  vc)  Hugoniot.  In  the  limit,  if  the  specific  volume  v  is 
reached  through  an  infinite  number  of  infinitesimal  shocks,  it  should  not  be  surprising  that  the  path  of 
compression  will  lie  along  the  iscntropc  through  the  original  (p0,  vj  point,  since  an  infinitesimal  shoJc 
was  shown  previously  to  be  iscntropic. 

In  practice,  only  a  single  Hugoniot  is  experimentally  determined  to  characterize  the  behavior  of  a 
material— that  Hugoniot  which  has  as  its  (p„,  vj  state  ambient  conditions  of  pressure  and  specific 
volume.  Thus,  references  to  "the"  Hugoniot  for  a  material  refer  to  a  particular  ’’reference”  Hugoniot 
which  has  as  its  origin  the  ambient  (p„,  vc)  state. 


7.  EXAMINATION  OF  SEVERAL  HUGONIOT  FORMS 

7. 1  Constant  Sound  Sneed.  Constant  sound  speed  implies  that  all  disturbances  travel  at  the  same 
velocity,  or 
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Determining  the  particle  velocity  from  continuity  gives 


up  =  Ca(J.  pjp)  . 


It  is  clear  that  the  particle  velocity  can  never  exceed  the  shock  velocity  and  simultaneously  satisfy 
continuity,  so  the  implication  is  that  the  maximum  possible  particle  velocity  is  also  the  shock  velocity 
C0,  and  this  only  occurs  as  the  shocked  density  becomes  infinite.  It  should  be  clear  from  momentum 
considerations,  therefore,  that  the  maximum  pressure  achievable  in  a  constant  sound  speed  medium  is 
p0  C02.  Indeed,  the  relation  is  given  by 


Pa  -  Po 


P  eCl» 
(1  +  P) 


7.2  The  Ideal  Gas  Form.  Much  of  the  work  done  in  thermodynamics  is  done  on  the  assumption 
that  the  working  medium  may  be  approximated  as  an  ideal  gas.  Unfortunately,  at  large  compressions, 
materials  do  not  generally  behave  in  an  ideal  fashion.  Nonetheless,  we  examine  the  Hugoniot  form  for 
the  sake  of  academic  interest. 

An  ideal  gas  is  one  which  obeys  the  EOS  p/p  =  RT,  where  R  is  the  gas  constant,  and  T  is  the  gas 
temperature.  Additionally,  an  ideal  gas  is  assumed  to  have  constant  specific  heats  cp  and  cv.  Under 
these  conditons  of  idealness,  E  =  T,  and  the  EOS  may  be  expressed  as  p/p  =  (R/cv)  E.  But 
R  =  cp  -  cv  by  definition,  and  the  ratio  Cp/cv  is  often  expressed  as  y,  so  that 

p  =  p  (y  - 1)  E  . 

Consider  the  shock  energy  equation,  and  eliminate  energy  terms  from  it  by  substituting  the  gamma  law 
from  above.  In  that  case. 


Pk 

P  U  - 1) 


Po 

P,(Y  -  1) 


(P„+Pa>  (1/p„-  1/P)  • 
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This  may  be  solved  for  ihe  Hugonioi  pressure  ph,  and  expressed  in  compression  form,  by  substituting 
(1  +  P)  for  p/p0,  to  give 

pjp  « 

h  °  (Y  -  1)  (Y  -  D(1  -  lO 


Expressed  in  a  (pj,  -  p0)  form,  it  becomes 


Pk-Po 


_ 

2  -  (Y  -  1)  pt  ' 


This  Hugoniot  form  is  seen  to  approach  infinite  pressure  when 


U  =  2  /  ( Y  -  1)  . 


This  limit,  expressed  in  terms  of  density,  is 


plpo 


YJLl 

T-l  ' 


7.3  The  Linear  Uf  -  t^  Form.  It  has  been  empirically  observed  for  many  materials,  over  a 
significant  range  of  pressures,  that  the  relationship  between  the  shock  and  particle  velocity  can  often 
be  fit  by  a  straight  line,  as  in 


U.  =  Cn  +  S  u. 


Thus, 


Ph  -  Po  *  Pc  (C,  +  Sup)  up  , 
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and  from  continuity,  particle  velocity  may  be  determined  as 


t)  C 

- - ;  .  where  *1  =  1-  PjP  • 

1  -  TJ  i 


Panicle  velocity  may  be  eliminated  from  these  equations  to  give 


Pk-Po 


P  pCU 

(1  -TlS)2  ' 


Expressing  this  in  compression  form  gives: 

p„c.!MW) 

P,~P° "  (i-(S-l),.)J  ’ 


It  can  be  seen,  for  S  >  1,  that  this  linear  Us  -  Up  Hugoniot  equation  approaches  infinite  pressure  at 
finite  compression  (p  =  1/[S  -  1]),  not  unlike  the  ideal  gas  EOS.  Otherwise,  for  S  <  1,  the  Hugoniot 
pressure  converges  to  pG  C02/(S  -  l)2  for  large  p.  Also,  the  special  ease  of  S  =  0  reduces  to  the 
previously  derived  situation  of  coastant  sound  speed.  This  finite  limiting  pressure  for  S  <  1  arises 
because  the  U,  function  intersects  the  up  function  at  Us  =  CJ{\  -  S).  Since  U,  must  always  exceed  up. 
these  limiting  values  of  U$  and  up  cause  the  limiting  value  in  pressure  loo.  Though  this  U5  =  up 
condition  is  unrealistic,  apparently  several  materials  exhibit  S  <  1  at  lower  compressions  (Kohn  1969). 
Such  data  can  be  misleading,  and  can  cause  problems  if  S  <  1  is  retained  out  to  larger  compressions  In 
numerical  simulations. 


7.4  The  Polynomial  Hugoniot.  Often,  a  form  is  not  assumed  on  the  U,  -  Up  relation,  or  the 
compressibility  at  all,  but  rather,  an  empirical  fit  is  made  directly  to  the  Hugoniot,  as  in 

Pk'Po  =  £  Ki  P'  » 

where  Kj  are  the  empirical  constants,  and  n  is  the  order  of  the  polynomial  fit.  A  popular  value 
for  n  is  three,  giving  a  cubic  fit.  Note  the  functionally  different  forms  between  this  Hugoniot  form 
and  the  linear  U,  -  up  form,  for  instance.  To  work  the  problem  in  reverse,  and  obtain  the  shock 
velocity  for  a  cubic  Hugoniot  fit,  for  example,  one  can  manipulate  continuity  and  momentum 
cquatioas  to  get: 
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v)  =  [Kx  ♦  {Kx  +  KJ  n  ♦  (K*  +  K,  )  p2  ♦  K^)  /p.  . 


Note  that  Us:  at  p.  =  0.  which  is  CQ2,  equals  K,/p0. 


8.  THE  IMPORTANCE  OF  DISTURBANCE  VELOCITY  ON  NUMERICAL  STABILITY 

The  primary'  emphasis  in  describing  the  shock  loading  behavior  of  materials  has  been  on  the 
Hugoniot  fit  itself,  the  ph  vs.  p  relation.  However,  the  importance  of  the  disturbance  speeds  U, 
and  C  (i.c..  the  shock  velocity  and  local  sound  speed)  cannot  be  underestimated.  In  impact 
computations,  where  the  wave  motions  arc  the  important  part  of  the  problem  solution,  knowledge  of 
disturbance  velocities  becomes  necessary  for  accuracy,  and  in  explicit  schemes,  for  stability,  of  the 
calculation. 

In  explicit  integration  schemes,  which  arc  by  far  the  most  prevalent  in  use  for  impact  calculations, 
stability  requirements  place  limitations  on  the  integration  timestep.  In  particular,  the  explicit 
integration  timestep  must  be  limited  so  that  a  disturbance  can  traverse  no  more  than  a  single 
computational  zone  per  integration  cycle.  To  enforce  this  requirement,  one  needs  knowledge  of  both 
the  particular  zone  size  in  question  and  the  disturbance  velocity  in  that  computational  zone.  The 
former  is  simply  a  function  of  zone  geometry,  while  the  latter  is  the  shock  velocity  Uf,  when  the 
material  is  compressing,  and  the  sound  speed  C,  when  the  material  is  expanding.  Thus,  an  accurate 
knowledge  of  disturbance  velocities  U,  and  C,  under  various  states  of  loading,  is  vital  to  the  numerical 
stability  of  explicit  integration  schemes. 

In  addition  to  the  effect  of  the  disturbance  velocities  on  the  stability  of  the  explicit  numerical 
integration  scheme,  these  velocities  arc  intimately  related  to  the  stability  of  the  EOS  calculation  itself, 
even  in  the  absence  of  a  numerical  code  application.  There  arc  thermodynamic  considerations  which 
must  be  satisfied  by  any  valid  EOS. 

Three  particular  modes  of  EOS  instability  will  be  addressed  in  the  following  sections.  Stated  in 
English,  these  stability  criteria  sound  almost  trivial,  since  they  immediately  result  from  fundamental 
thermodynamic  statements,  which  have  been  known  and  professed  for  years.  Nonctliclcss,  these 
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criteria,  which  can  be  expressed  mathematically,  and  are  tested  in  the  framework  of  the  Mie-Grilneisen 
EOS,  show  a  necessary  interrelation  amongst  various  EOS  parameters,  which  is  not  commonly 
appreciated.  As  such,  indiscriminate  choice  of  parameters  to  drive  the  EOS  will  often  cause  at  least 
one  of  these  criteria  to  be  violated.  Use  of  such  data  in  a  numerical  scheme  can  result  in  bizzarre 
manifestations  of  numerical  ins'ability. 

The  three  stability  criteria  are  given  below  and  apply  to  materials  where  phase  changes  and/or 
porosity  are  not  considered  (i.c.,  simple  substances).  For  the  first  two  criteria,  consider  a  state  (p,  v), 
which  lies  on  the  Hugoniot  originating  from  the  state  (p0,  v0). 

(I)  at  the  state  (p,  v),  the  slope  of  the  Hugoniot  with  respect  to  density  must  exceed  the  slope  of 
the  Rayleigh  line  connecting  states  (pQ1  v0)  and  (p,  v); 

(II)  at  the  state  (p,  v),  the  slope  of  the  isentrope  must  lie  between  the  slope  of  the  Hugoniot,  and 
the  slope  of  the  Rayleigh  line  connecting  states  (p0,  v0)  and  (p,  v); 

(III)  the  slope  of  an  isentrope  with  respect  to  density  is  always  greater  than  zero. 

From  the  second  law  of  thermodynamics,  it  is  known  that  entropy  is  monotonically  increasing  along  a 
Hugoniot.  By  definition,  entropy  is  constant  along  an  isentrope.  Finally,  a  Rayleigh  line  is  known  to 
have  one  point  of  maximum  entropy.  It  can  be  shown  (Zuker  1977)  that  this  point  of  maximum 
entropy  lies  between  the  states  (p0,  v0)  and  (p,  v).  Thus,  at  the  state  (p,  v)  and  traveling  in  the 
direction  of  increasing  density,  entropy  must  be  decreasing.  Criteria  I  and  II  express  this  relational 
placement  of  the  Rayleigh  line,  the  isentrope  and  the  Hugoniot.  Criterion  III  simply  implies  the 
obvious,  that  the  speed  of  sound  in  a  material  is  a  positive  quantity. 


9.  MODE  I  INSTABILITY  OF  SEVERAL  HUGONIOT  FORMS 

It  was  mentioned  that  a  Hugoniot  has  the  slope  (with  respect  to  density  p  or  compression  p)  of  the 
isentrope  only  at  me  end  point  of  the  Hugoniot.  At  other  points  along  the  Hugoniot,  the  Hugoniot  will 
always  be  steeper  than  the  corresponding  isentrope.  Let  us  examine  the  two  most  popular  Hugoniot 
forms  (linear  U,  -  up  fit  and  cubic  fit)  to  check  their  adherence  to  these  rules.  The  slope  of  the 
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Hugoniot  at  any  point  on  the  Hugoniot  is  simply  dp^/dp.  Total  derivatives  are  used  since,  along  the 
Hugoniot,  pressure  is  only  a  function  of  a  single  variable  (p,  p,  or  v,  at  the  discretion  of  the  user).  To 
determine  the  slope  of  the  isentrope  at  that  point,  one  can  make  use  of  the  definition  that  dp/dp  I , 
equals  the  square  of  the  sonic  velocity. 

For  the  case  of  the  linear  Us  -  up  lit,  recall  that 

Ph'P°  (1  -(5-  l)p)2  ' 


Taking  the  derivative  with  respect  to  p,  one  may  obtain 

dp,  _  P „cl  [1  +  (S  +  i)n] 

~dP  [1  -  (S  -  1) 


The  chain  rule  dictates  that  dp^dp  =  dp^dp  dp/dp.  where  dp/dp  is  simply  l/p0,  so  that 

_  c,2  [i  +[s  +  i)p] 

dp  "  [1  -  (S  -  1)  p]J 


Without  knowing  additional  information  about  the  F.OS,  the  slope  of  the  post  shock  isentrope 
3p/dp  I  s  =  C2,  and,  thus,  the  sonic  velocity  C,  cannot  be  explicitly  determined  after  a  shock. 

However,  the  particle  velocity  relative  to  the  stationary  shock  front  after  a  shock,  equal  to  (U,  -  Up),  is 
known  to  be  subsonic.  It  may  also  be  shown  that  dp/<?p  along  the  Rayleigh  line  connecting  the  pre- 
and  post-shocked  states  exactly  equals  (U,  -  u^)2.  If  the  Hugoniot  slope  is  larger  than  the  slope  of  the 
actual  isentrope,  then  clearly  its  slope  must  be  larger  than  that  of  the  Rayleigh  line,  which  represents 
the  lower  bound  on  the  slope  of  the  isentrope. 

From  continuity,  the  quantity  (U,  -  Up)  may  be  expressed  either  as  U/(J  +  p)  or  Up/p.  Thus,  ihe 
slope  of  the  isentrope  dp/dp  1 1  =  C2  mu  it  be  greater  th?n  U,2/(l  +  ft)2,  which  is  the  slope  of  the 
Rayleigh  line.  Now,  consider  the  Rayleigh  line  slope  at  some  point  on  the  Hugoniot.  For  the  linear 
Us  -  Up  assumption,  Us  =  CQ  +  S  up.  However,  it  is  desired  to  obtain  U,  in  terms  of  compression  p. 
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Elimination  of  up  from  the  continuity  equation  gives. 


P  = 


U  -u 
*  p 


U,-(U,-C0)IS  ' 


Upon  solution  for  U;, 


U . 


CJl  +  y i) 


Thus.  dp/9p  1R  =  Us2/(1  +  p)2.  and  is  given  as 


dp  _  _ _ 

ap  je "  n-(5-Dp]2 ' 


For  the  Hugoniot  slope  to  exceed  the  corresponding  Rayleigh  line  slope  for  all  p  >  0.  the  ratio  of 
dp^dp  to  3p/3p  I R  must  always  exceed  unity.  Failure  to  do  so  implies  that  there  can  exist  no  valid 
equation  of  state  to  describe  the  material.  We  will  call  Hugoniots  whose  slope  docs  not  exceed  that  of 
the  corresponding  Rayleigh  line,  as  Mode  I  unstable.  Mathematically,  the  condition  is  stated  below: 

MODE  I  CRITERION: 


dpjdp  ^ 
dp/d  pj. 


1 ,  for  p  >  0,  along  the  Hugoniot. 


For  the  linear  U,  -  up  form,  the  stability  ratio,  which  must  exceed  unity,  is  given  by 


dpjdp  _  [l  +(S+  l)p] 
dpidp  |,  "  [1  -  (5  -  1)  p] 


>  i ,  for  p  >  0,  along  the  Hugoniot. 


Clearly,  for  any  S  greater  than  zero,  the  stability  criterion  is  not  violated.  Thus,  the  linear  U,  -  Up 
formulation  is  always  Mode  I  stable. 
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Considering  now  the  cubic  empirical  Hugoniot  fit,  it  is  of  interest  to  see  how  this  Hugoniot  form 
fares  with  the  Mode  I  stability  criterion.  Mimicking  the  operations  performed  on  the  U,  -  Up  forms, 
dp^dp  is  determined  as 

£.[*,.  2*  P.. 

The  slope  of  the  Rayleigh  line,  dp/9p  I R  equals  Us2/(1  +  p)2,  and  is  given  by 

Kt  *  p  +  (*2  +  *3)  *  *3  ^ 

P0  O  + 

Finally,  the  Mode  1  stability  criterion,  which  insures  that  the  slope  of  the  Hugoniot  exceeds  that  of  the 
Rayleigh  line  connecting  the  pre-  and  post-shocked  states,  is  given  as 

dpjdp  _  (K,*2 K,  p  +  3K,  p2)(l  +  p)2  >1 

dp/d  pi,  '  Kx  *  (Kx  +  K2)vl  +  +KJ  p2  +  K,  p3 

for  p  >  0,  along  the  Hugoniot.  Because  of  the  dependence  of  this  term  on  three  independent 
parameters,  a  simple  answer  cannot  be  given.  However,  the  data  collected  and  fit  by  Kohn  (1969)  into 
K,.  K2,  K3  form  can  still  be  checked  for  slope  violations.  It  was  found  that  21  of  77,  or  27%,  of  the 
materials  violated  the  criterion  in  part  of  the  compression  range  where  the  empirical  fit  was  supposedly 
valid.  Perhaps  this  behavior  can  be  attributed  to  the  forcing  of  a  cubic  Hugoniot  function  onto  a 
highly  non-linear  material  response,  possibly  due  to  phase  changes.  In  such  cases.  Mode  I  instability 
implies  that  a  cubic  Hugoniot  fit  is  inadequate.  The  ramifications  of  Mode  I  instability  on 
computational  stability  arc  no:  as  clear  as  the  other  modes  to  be  discussed,  where  instability  is  more 
likely  to  guarantee  computational  catastrophe. 


0£ 

dp 


10.  INTRODUCTION  TO  THE  MIE-GR0NEISEN  EQUATION  OF  STATE 

Traditionally,  an  EOS  is  an  equation  which  describes  the  pressure  of  a  material  as  a  function  of 
density  and  temperature.  It  is  needed  to  solve  the  continuity,  momentum,  and  energy  equations  which 
govern  the  thermodynamic  transition  of  a  material.  In  the  case  of  impact  though,  thermodynamic 
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states  change  so  rapidly  that  there  is  little  time  for  heat  transfer.  Under  these  conditions  of  adiabatic 
transition,  the  energy  equation  need  not  explicitly  contain  temperature,  and  the  governing  equations 
may  tie  solved  if  cither  the  entropy  or  the  internal  energy  is  known  as  a  function  of  pressure  and 
dcasity  (Zcldovich  and  Razier  1966).  The  Mic-Griinciscn  EOS  is  of  this  variety,  in  that  internal 
energy,  and  not  temperature,  is  related  to  pressure  and  density.  The  penalty  for  not  having  the 
traditional  p  (p,  T)  form  lation  is  that  temperature  is  not  a  computable  state  variable. 

For  impact  calculations,  accurate  description  of  shock  transition  states  is  very  important. 
Additionally,  high  pressure  material  data  arc  usually  obtained  from  shock  transition  experiments  and  is 
given  in  the  form  of  a  Hugoniot  curve.  The  Hugoniot  provides  the  means  to  ascertain  the  pressure 
and  internal  energy  as  a  function  of  compression  for  those  thermodynamic  states  which  can  be 
achieved  through  shock  compression  of  a  material. 

Thus,  realizing  that  an  EOS  is  only  an  approximation  to  the  actual  behavior  of  a  material,  it  seems 
desirable  to  have  an  EOS  which  is  very  accurate  along  the  path  where  data  arc  available  (i.c.,  along 
the  experimentally  derived  Hugoniot).  State  points  other  than  those  found  directly  on  the  Hugoniot  are 
obtained  from  an  EOS  which  references  material  back  to  the  state  along  this  reference  Hugoniot.  The 
Mic-Griineisen  EOS  provides  the  framework  through  which  state  points  can  be  tied  back  to  an 
arbitrary  reference  function,  as  described  below. 

Consider  the  textbook  form  of  the  Mic-Grilnciscn  EOS  (McQuarric  1976),  expressed  in  terms  of 
dcasity 

P  =  P2  ~  +  rP(£-tf>  ♦ 

dp 

where  U  is  die  specific  interatomic  potential  energy  (a  function  of  density  only),  and  E  is  the  specific 
internal  energy  comprised  of  both  interatomic  potential  (i.c.,  compressive)  energy  and  atomic 
vibrational  (i.c.,  thermal)  energy.  Note  that  the  term  p2  dU/dp  represents  the  "cold  pressure."  which  is 
that  pressure  arising  from  compression  at  absolute  zero  temperature,  in  the  absence  of  specific 
vibrational  internal  energy  (E-U).  The  term  T  is  the  Grilnciscn  coefficient,  which  is  a  function  of 
density  only,  derived  by  assuming  a  particular  form  on  the  vibrational  frequencies  of  the  crystal  lattice. 

If  the  cold  pressure  curve  is  known  (or  assumed),  and  if  the  thermodynamic  states  arc  known 
along  another  reference  curve  fe.g.,  the  Hugoniot),  then  the  Ortlnciscn  coefficient  may  be  inferred 
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directly,  as  a  function  of  density.  On  the  other  hand,  if  the  thermodynamic  states  arc  known  only 
along  a  "non-cold"  reference  curve,  but  if  a  functional  relationship  is  assumed  for  the  Griinciscn 
coefficient,  then  the  interatomic  potential  function  may  be  implicitly  removed,  by  referencing  the 
thermodynamic  state  at  any  point  (p.  p,  E)  back  to  the  state  along  the  reference  curve,  at  the  same 
density  (p.c  f  P*  Efcf)’  This  is  done  by  taking  the  difference  of  the  Mic-Griinciscn  EOS,  given  above, 
at  the  two  state  points  in  question.  It  is  this  approach  which  is  almost  universally  adopted  by 
hydrocodcs,  to  give 

(p  -  prfJ)  =  r  p  (E  -  . 

It  should  be  noted  that  die  interatomic  potential,  U,  and  thus  the  cold  pressure  curve,  may  be 
recovered,  by  solving  the  following  non-linear  ordinary  differential  equation: 

du  r u  _  <t>(p) 

dp  P  p2 


where  0(  p)  is  a  function  of  density  only,  obtained  from  the  "non-cold”  reference  curve  and  the 
assumed  Griinciscn  function: 

4>(p)=p^-rp£^. 


For  the  present  purposes,  the  reference  function  is  chosen  as  die  Hugoniot,  so  that  the  Mic-Griinciscn 
EOS  becomes 

(P-Ph)  =  Tp(E  -EJ  . 

The  quantities  Pj,  and  Eh  arc  the  pressure  and  specific  internal  energy  of  the  reference  Hugoniot  state 
point  at  the  same  density  as  the  state  (p,  p,  E)  in  question.  Often,  the  shock  energy  equation  is 
substituted  directly  into  the  Mic-Griinciscn  EOS  to  give 

P  =  />*(!  -Fii/2)  ♦  p,T  (1  ♦  y)  [E-E0]  -  p0Ypl2  . 


The  ambient  conditions  p0,  Ec,  arc  generally  pegged  at  zero,  so  that  the  Mic-Griinciscn  EOS,  with 
Hugoniot  reference,  is  expressed  as 

p  <1 -I>/2)  ♦  p,r(l4p)£. 
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Note,  in  an  analogous  fashion  to  the  Mic-Grtinciscn  EOS.  that  ideal  gas  behavior  may  be  modeled 


,  given  by  f  for  solids,  is  held  constant  for  gases,  equal  to  (y  -  1),  and  if  the  Hugoniol 


fonm  is  that  for  an  ideal  gas;  namely. 


Pk  ~Po 


2yn 

2-(y-i)n  * 


The  Mic-Griincsicn  EOS,  with  Hugoniol  reference,  is  only  valid  where  the  Hugoniot  reference  is 
valid.  If  the  reference  Hugoniol  is  obtained  by  shocking  material  from  ambient  conditions,  then 
clearly  the  Hugoniot  is  only  valid  in  some  region  p  >  0.  Hydrocodc  implementors  often  fail  to  deal 
with  this  fact  accurately.  It  is  well  known  that  the  Hugoniot  curve  is  a  very  close  approximation  to 
the  isentrope  at  small  values  of  compression,  and  so  the  pressure  reference  function  usually  employed 
by  hydrocodcs  in  tension  is  the  isentrope,  which  is  generally  approximated  by  a  curve  which  matches 
value  and  slope  of  the  Hugoniot  at  p  =  0.  There  is  nothing  wrong  with  this  assumption,  but 
innaccuracics  arc  introduced  when  the  shock  energy  equation  is  still  employed  directly  to  compute  the 
energy  reference  function,  as  in  EPIC  (Johnson,  Colby,  and  Vavrick  1978),  which  is  tantamount  to 
assuming  the  existence  of  a  tensile  shock. 

Clearly,  the  energy  reference  function  to  drive  the  Mic-Griinciscn  EOS  in  tension  should  not  be 
the  shock  energy  equation 


E„  -  Eo  - 


P*P 


2p.(l  +  l x) 


(where  p  ~  0)  , 


but  rather  some  physically  plausible  function — for  instance,  the  iscntropic  expansion  integral  of  work: 


Enf~  Ec~  /  Pnf 


Pc  (1  +  P)2 


One  possible  approximation  for  the  iscntropic  release  function  for  p  <  0  is  to  assume  constant 
compressibility.  In  this  ease,  where  pc  =  E0  =  0  is  assumed  for  simplicity. 
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JW  =  KeV  *  and 


£*  =  [p« - =  Wp.)  • 

O  PcH-n)2 


Comparing  ihc  internal  energy  from  this  reference  function  to  that  computed  from  the  invalid 
application  of  the  shock  energy  equation  in  the  tensile  region,  given  by 

E  -  -  K°  ^ 

A  2  pfl ( 1  +  p)  2p#(h(i)  ’ 

one  may  observe  that  the  difference  in  computed  internal  energies  is  significant,  as  shown  in  Table  1. 


TABLE  1.  Comparing  Internal  Energies  Resulting  From  Constant 

Compressibility  "Iscntropic"  ar.d  "Tensile  Shock"  Expansions 


T - 

It 

v/v„ 

Po  E/K0 
(iscntropic) 

Po  E|/K0 
(tensile  shock) 

%  Error 

r  -■  ■  - ‘-1-  ■ - — 

().() 

1.00 

0. 

0. 

0. 

_  2 

1.25 

0.027 

0.025 

-7.4 

-.4 

1.67 

0.156 

0.133 

-14.7 

-.6 

2.50 

0.584 

-23.1 

-.8 

5.00 

2.39 

1.60 

-33.1 

Oilier  reasonable  choices  of  iscntropic  release  pressure  reference  functions  may  also  be  used  in  the 
region  p  <  0.  For  example,  let  the  reference  sound  speed  for  a  material  in  tension  (p  <  0)  be 
described  with  the  following  function: 


For  this  particular  sound  speed  reference  function,  the  tangent  bulk  modulus  may  be  computed  as 

K  =  Ke{\ ♦  p)^°>  . 
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This  function  has  many  nice  properties.  Two  such  advantages  are  that  it  reduces  to  constant 
compressibility  for  g  =  -1/2,  and  to  constant  sound  speed  for  g  =  0.  Additionally,  the  sonic  velocity  is 
at  the  proper  value  at  p  =  0  and,  for  g  greater  than  zero,  the  sonic  velocity  drops  off  to  zero  when  the 
material  is  infinitely  expanded  (at  p  =  -1).  Increasing  the  exponent  g  causes  the  sonic  velocity  to  drop 
off  more  rapidly  in  tension  and  will  limit  the  maximum  value  of  tension  achievable  during  an 
iscntropic  expansion,  which  is  a  material  feature  observed  experimentally  and  often  modeled 
computationally  by  way  of  a  pmin  parameter.  The  reference  pressure  function  is  acquired  by 
integrating  the  tangent  bulk  modulus  and  is  given  below  for  g  not  equal  to  (-1/2): 


P 

p«-!* 

o 


dy. 

(1  +  p)2 


K. 

(2  g  *  2) 


[(1  -  1]  . 


The  reference  energy  function  (for  g  not  equal  to  -1/2)  is  acquired  by  integrating  the  reference  pressure 
function: 


'r* 


=  Jp« 


d\i 


P.O  +  P)3 


-  K. 


[(1  +  p)(2«‘2)  -  1]  -  p (2#  +  2) 

P,(2g  +  2)  (2g  +  1)  (1  +  n) 


For  the  special  case  of  constant  sound  speed  adiabatic  expansion  from  p  =  0,  the  exponent  g  equals  0, 
and  one  gets: 


Pn,  -  *„(p  +  p2)/ 2  ,  and 


Kg  p2 

2  pe  (i_4  p) 


On  the  other  hand,  the  reference  pressure  may  be  limited  to  a  tensile  pressure  of  pmij),  which  occurs  at 
infinite  expansion,  by  defining  the  exponent  g  to  be 


8  = 


2  IPtninl 


-1  . 
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Using  steel,  as  an  example,  where  Kc  =  1,380  kbar,  and  pmin  =  -35  kbar,  g  takes  on  a  value  in  the 
vicinity  of  19.  For  aluminum,  where  K0  -  700  kbar,  and  pmin  =  -10  kbar,  the  value  g  takes  is 
approximately  34.  Table  2  depicts  the  parameter  relationships  for  g  =  20. 


TABLE  2.  Values  of  Compression.  Relative  Volume,  and  Normalized  Values  of  Sonic 

Velocity,  Pressure,  Tangent  Bulk  Modulus,  and  Secant  Bulk  Modulus,  Respectively, 
for  the  Proposed  Tension  Limiting  Model,  With  g  =  20 


v/v0 

c/c0 

Pre/Pmin 

K/K0 

P/K0P 

0.00 

1.00 

1.000 

0.000 

1.000 

N/A 

-.002 

1.002 

.961 

.081 

.921 

.960 

-.02 

1.02 

.668 

.572 

.437 

.681 

-.04 

1.04 

.442 

.820 

.188 

.488 

-.06 

1.06 

.290 

.926 

.079 

.367 

-.08 

1.19 

.189 

.970 

.033 

.289 

-.10 

1.11 

.122 

.988 

.013 

.235 

1 1.  THE  GRUNEISEN  PARAMETER 

The  Grunciscn  parameter  T  is  an  important  parameter  used  by  the  Mic-Griincisen  EOS.  Though 
the  assumption  that  T  is  a  function  of  volume  only  is  only  an  approximation,  the  parameter  itself  is 
not  empirical  but  in  fact  carries  physical  significance.  For  ideal  gases,  the  parameter  analogous  to  T 
has  been  shown  to  be  (y  -  !)•  For  solids,  the  parameter  also  has  meaning. 

To  compute  its  value,  consider  a  cube  of  material  at  state  1  characterized  by  pressure  p  and 
volume  v.  Let  the  material  undergo  a  two-stage  thermodynamic  process.  First,  add  heat  dQ  at 
coastant  pressure  in  order  to  go  to  state  2,  characterized  by  p,  V  +  dV.  Finally,  and  with  no  heat 
addition,  isentropically  compress  the  material  back  to  its  original  volume  in  order  to  arrive  at  state  3, 
characterized  by  p  +  dp,  V. 

During  the  process  1-2,  the  cube  of  material,  the  length  of  whose  side  was  originally  L,  expands 
so  that  the  length  is  now  L  (1  +  adT),  where  a  is  the  coefficient  of  linear  expansion.  Thus,  the 
relative  volume  change  dV/V  is  given  by  3adT.  The  amount  of  heat  added  per  unit  mass  dQ.  is 
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cp  dT.  For  solids,  which  are  nearly  incompressible  in  comparison  to  gases,  the  specific  heats  at 
constant  pressure  and  volume,  cpl  and  cv,  arc  nearly  identical.  Thus,  the  heal  added  can  be 
approximated  by  cv  dT.  The  temperature  change  dT  may  be  eliminated  by  relating  dT  back  to  volume 
change  dV  so  that 


dQ 


cv  dV 
3o  V  ' 


The  work  done  by  the  material  per  unit  mass  is  negative  since  the  volume  change  is  the  opposite  sense 
of  the  pressure  loads  and  is  given  by  dW  =  -p  dV/(pV).  From  the  first  law  of  thermodynamics,  the 
internal  energy  change  dE,  given  by  dQ  -  dW.  can  be  expressed  as 


pdV 
pV  ' 


The  process  from  state  2  to  3  has  no  heat  addition,  so  dQ  =  0.  The  work  done  by  the  system  is 
dW  =  (p  +  dp/2)  dV/(pV)  which,  to  the  first  order,  is  p  dV/(pV).  The  rise  in  internal  energy  is 
therefore  the  negative  of  the  work 


p(-dV) 

pV 


Thus,  the  difference  between  the  pressures  going  from  slate  1  to  3  is  dp,  and  the  rise  in  internal 
energy  is 


dE 


13  ~ 


crdV 
3a  V  ' 


The  volume  of  states  1  and  3  are  identical,  so  that  the  Grvinciscn  parameter  T  will  remain  constant  for 
the  direct  transition  from  state  1  to  3.  when  using  the  Mic-Griinciscn  EOS,  which  relates  states  1  and  3 
as  follows: 


But  path  23  is  iscntropic,  so  that  dp/C-dV)  is  dp/dV  |  s  =  (-p/V)  Op/dp  I Bui  <)p/3p  I ,  is  simply  the 
square  of  the  sound  speed  C,  or  alternately,  the  ratio  of  the  bulk  modulus  and  the  density  K/p. 
Substuting,  and  solving  for  T  gives 


3a  C2  =  3a  K 
Cv  P  Cv 


From  the  Mic-Griinciscn  EOS,  it  is  clear  that  the  Griinciscn  parameter  T  relates  changes  in  internal 
energy  to  changes  in  pressure  at  constant  volume.  As  can  be  seen  from  the  dimensional  analysis  of 
the  second  expression  for  T  which  follows,  the  numerator  expresses  the  pressure  rise  per  increment  of 
temperature,  while  the  denominator  expresses  the  (volume-based)  energy  rise  per  increment  of 
temperature.  These  terms  combine  in  a  way  to  produce  a  non-dimensional  expression  which 
represents  the  Griinciscn  parameter. 

dVjV  dp 

r  =  Aa_*  =  deg  dVjV  =  dp!  deg  =  f 

P  cv  dE_  p  dEldeg  pdE 

deg 

Note:  As  in  oilier  places  in  this  document,  dE  represents  energy  change  per  unit  mass. 


12.  MODE  II  INSTABILITY  OF  THE  MIE-GRUNEISEN  EQUATION  OF  STATE 

Though  there  is  nothing  inherently  wrong  with  the  Mic-Griinciscn  EOS,  there  arc  often  situations 
where  inappropriate  choice  of  data  forces  a  thermodynamic  violation.  As  pointed  out  in  the  previous 
section  dealing  with  Mode  I  instability,  there  arc  certain  empirical  Hugoniot  fits  which  may  violate 
fundamental  thermodynamic  rules.  Dearly,  plugging  such  fits  into  Mic-Griinciscn  makes  the 
possibility  of  catastrophic  EOS  failure  significant.  Additionally,  there  arc  inconsistent  combinations  of 
Hugoniot  and  Griinciscn  functions  which  cause  thermodynamic  violatioas  which  are  not  immediately 
apparent. 

To  study  these,  consider  a  material  in  the  state  p,  Eh  on  the  Hugoniot  originating  at  (p0 ,  v^. 
From  this  state,  consider  two  possible  paths  resulting  from  a  small  compressive  increment,  dp:  the 
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path  along  the  isentrope  and  also  the  path  along  the  Hugoniot.  Along  the  isentropic  path,  the  pressure 
change  is 


dp 

dp 


dp  -  C2  d p  . 


The  internal  energy  change  along  the  same  path  is  the  negative  of  the  work,  or 


a e 

dp 


dp  -  Ph  (dplp2)  . 


Along  the  Hugoniot,  the  pressure  change  is  simply  (dp^/dp)  dp,  and  the  energy  change,  may  be 
obtained  by  taking  the  derivative  of  the  shock  energy  relation  -  E0  =  1/2  (p,,  +  pc)  (l/p0  -  1/p), 
to  get 


dp  =  1/2  [(ph+pe)lpl  +  (1/P„-  1/P)  dpjdp]dp  . 


If  this  material  is  to  obey  the  Mic-Grunciscn  EOS,  then  at  the  density  p  +  dp. 


„  dp 

\  op 


dp 


dph  , 

ph  ♦  --  efp 


=  r  p 


l  *  dp 


dp 


£.  +  - 


dE. 


dp 


dp 


Canceling  and  substituting  terms  gives 

(c2  ■  dPh/dp)  =  r  p  (p./p2  - 1/2  i(Ph  +  pjip2  +  (i/p0  -  up)  dPh/dp]) . 

Expressing  this  in  terms  of  compression  p,  and  solving  for  T  gives 

r  -  -  ldPk!dp  -  pp0c2] 

P  [p  dpjdp  -  (ph  -Pe)l(l  *  p)] 
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The  significance  of  this  expression  is  that,  for  a  given  Hugoniol  relationship  Pb  as  a  function  of  p,  the 
Griinciscn  coefficient  P  and  die  sound  speed  C  arc  related.  Since  the  sound  speed  is  known  to  have 
thermodynamic  limitations  placed  upon  it,  these  limitations  translate  into  constraints  on  the  Griineisen 
parameter  P.  Use  of  the  Mic-Griincisen  EOS  in  codes,  without  the  appropriate  restraints  on  T,  will 
usually  result  in  catastrophic  failure  of  die  EOS. 

To  examine  these  constraints  on  P.  consider  the  maximum  and  minimum  possible  values  on  C. 
Recall  dial  die  post  shock  particle  velocity  relative  to  the  shock  front,  equal  to  (Uf  -  Up),  is  known  to 
be  subsonic,  but  approaches  sonic  conditions  for  infinitesimal  shocks.  Thus,  as  an  absolute  lower 
bound,  let  C  approach  (U,  -  Up).  This  choice  of  C  gives  an  upper  bound  on  T.  But  (U,  -  Up)  equals 
Us/(1  +  p),  so  that  (U,  -  Up)2  is  given  by 

(£/,-«„)2  =  u)n  i  +  p)2  =  (ph-p0)i[po\i(i  +  p)i . 


Substituting  this  into  die  expression  for  T  gives,  as  an  upper  bound, 

r  =  2/p  . 


For  the  lower  bound  on  T,  one  may  assume  that  the  isentrope  approaches  the  HugonioL  In  this  case, 
dp/dp  |  s  =  clpf/dp.  Since  dp/<)p  I ,  equals  p„  C2  by  definition,  the  minimum  value  on  T  may  be 
determined  as  zero.  By  combining  these  two  criteria,  and  calling  them  the  Mode  II  stability  criterion, 
one  gets 

MODE  II  CRITERION: 


0  <  T  <  (2/p)  . 

Failure  to  satisfy  this  criterion  implies  that  the  slope  of  the  post  shock  isentrope  docs  not  fall 
between  that  of  the  Rayleigh’  line  and  the  Hugoniol. 

Recall,  from  the  ideal  gas  form,  that  the  parameter  analogous  to  T  is  constant,  and  equal  to 
(Y  -  1).  it  may  at  first  appear  that  the  ideal  gas  form  violates  (he  Mode  II  stability  criterion  for  laigc  p. 
However,  recall  that  infinite  pressure  was  approached  when  the  ideal  gas  compression  p  equaled 
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2/(y  -  1).  which  is  2 IV  in  Mic-Griineiscn  terms.  Thus,  over  the  physical  p  domain  of  the  ideal  gas 
form,  which  is  p  <  2/T,  Mode  II  stability  is  not  violated. 

In  general  though,  it  becomes  clear  that  the  terminal  value  of  the  Griineisen  coefficient,  let  us  call 
it  rx.  must  be  less  than  or  equal  to  2/px,  where  px  is  defined  as  that  value  of  compression  which 
produces  infinite  pressure  in  the  Hugoniot.  For  the  linear  Us  -  up  form,  px  occurs  when  the 
denominator  (1  -  (S  -  1)  p)  becomes  zero,  which  occurs  only  if  S  exceeds  unity.  In  this  ease, 


=  1/(S-1). 


The  value  of  rx,  therefore,  should  be  less  than  or  equal  to  2(S  -  1).  For  S  <  1,  the  linear  Us  -  u(> 
form  becomes  unbounded  only  as  p  approaches  the  infinite.  In  this  ease,  Tx  must  approach  zero. 
Similarly,  for  polynomial  fit  Hugoniots,  shock  pressure  becomes  unbounded  only  when  compression 
approaches  infinity.  Thus,  the  Griineisen  coefficient  T  should  also  approach  zero,  all  die  while 
remaining  less  than  2/p,  as  compressions  become  very  large. 

It  is  unfortunate  that  data  for  F  arc  generally  only  available  at  the  ambient  p  =  0  condition.  As  a 
result  of  this,  it  has  been  common  practice  to  leave  the  value  of  T  constant  in  numerical  computations 
for  lack  of  better  data.  In  analyzing  the  data  of  Kohn  (1969),  it  is  observed  that  none  of  the  data 
presented  violate  the  Mode  II  stability  criterion,  0  <  T  <  2/p,  for  the  range  of  p  over  which  the  data 
arc  collected.  However,  should  the  data  be  extrapolated  to  larger  values  of  p,  Mode  II  instabilities 
may  result,  since  100%  of  the  cubic  fils,  and  74.5%  of  the  linear  U,  -  Up  fits  described  by  Kohn 
(1969)  will  violate  Mode  II  stability  at  larger  p  if  T  is  held  constant  at  its  initial  value. 


13.  CORRECTIONS  FOR  MODE  II  INSTABILITY 


To  enhance  Mode  II  stability,  some  code  developers  have  imposed  functional  forms  on  the 
Griineisen  coefficient.  In  the  case  of  the  EPIC  codes,  this  is  not  done,  and  T  is  held  constant  at  To. 
As  pointed  out.  Mode  II  criterion  will  be  violated  for  p  >  2/T0.  In  HULL  (Matuska  and  Osbom 
1987),  the  value  of  F  is  a  defined  as  T  =  ro  p^p,  which  in  tcims  of  compression  p,  is 


r  =  r0/(i  + ») . 


This  form  can  still  violate  the  Mode  II  criterion  if  ro  >  2  and  p  >  2/(F0  -  2).  DYNA  (Hallqutst 
1989),  as  well  as  versions  of  CALE  (Tipton  1989),  use  a  first  order  correction  to  T0,  which  is 
equivalent  to  assuming  the  form  on  the  Griineisen  parameter  as 


r  =  (r„+>4p)/(i  +  p) , 


where  A  is  a  parameter  which,  when  fit  to  EOS  data  at  moderate  pressures,  is  always  greater  than  or 
equal  to  zero  and  generally  lies  in  the  vicinity  of  0.5  for  most  metals  (Tipton  1989).  If  the  value  of  A 
is  chosen  as  zero,  the  DYNA/CALE  form  for  T  reduces  into  that  used  by  HULL.  Even  though  a 
positive  value  for  the  parameter  A  might  cause  the  data  to  fit  better  at  the  pressures  of  interest,  it  also 
has  the  side  effect  of  lessening  stability  at  larger  compression.  If  A  is  positive,  the  limits  on 
compression  which  assure  Mode  II  stability  arc 

(2-r,),|(2-ry,81f  A>0 

*  2  A 


For  materials  like  aluminum,  with  T0  =  2.09  and  A  =  0.49,  the  DYNA/CALE  Griineisen  form  has  a 
limiting  compression  of  stability  near  2.  This  value  is  better  than  the  limiting  value  of  unity,  derived 
from  the  constant  T  =  2  assumption,  but  is  still  prone  to  failure  when  simulating  hypcrvclocity  impact. 
MESA  (Bolstad  1990)  uses  a  form: 

r  «  Y«  ♦  Yi/Cl  ♦  lO  » 


which  is  directly  convcrtablc  into  the  CALE  version,  if  one  takes  ro  —  (y0  +  yt)  and  A  -  y0.  Thus, 
the  arguments  made  about  the  CALE  form  apply  to  MESA  as  well. 

From  the  perspective  of  assuring  Mode  II  stability,  even  at  very  large  compressions,  an  alternate 
Griineisen  form,  being  proposed  here  is  given  as  follows: 
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r  = 


rfl 

i  *  Pi*  * 


where  (3  is  the  newly  introduced  parameter.  Unlike  the  other  forms  on  F,  where  stability  is  a 
coincidental  by-product  of  the  choice  of  ro.  Mode  II  stability  can  always  be  ensured,  with  the  current 
form,  given  appropriate  choice  of  (3.  The  form  reduces  to  that  used  by  HULL  if  P  is  forced  to  unity 
and  can  be  made  to  approach  constant  T  if  p  equals  0.  However,  constraining  P  to  a  value  greater 
than  or  equal  to  YJ1  will  ensure,  for  any  Hugoniot  form,  that  Mode  II  stability  is  always  satisfied  for 
any  finite  p.  The  constraint  on  p  can  be  relaxed  somewhat  if  the  associated  Hugoniot  form  is  valid 
only  in  some  finite  range  0  <,  p  <  p,.  An  example  of  such  a  limited  domain  Hugoniot  is  the  U,  -  up 
form,  with  S  >  1,  where  p*  =  1/(S  -  1).  For  these  eases,  the  constraint  on  p,  to  ensure  Mode  II 
stability,  is 


P  * 


r,-(2/p,) 

2 


To  see  how  this  proposed  form  can  match  the  DYNA/CALE  form  at  lower  compressions  but  still 
provide  stability  at  higher  compressions.  Figure  2  is  provided.  In  it,  three  Griinciscn  relationships  for 
aluminum  are  shown  (data  for  aluminum,  T0  =  2.09,  and  S  =  1.33  arc  extracted  from  Kohn  [1969]). 
The  A  =  .49  (CALE)  curve  shows  the  functional  form  of  T,  using  the  DYNA/CALE  functional  form. 
The  value  of  A  =  .49  comes  directly  from  the  CALE  manual  (Tipton  1989)  and  was  fitted  from 
experimental  data.  At  large  compressions  (in  excess  of  2),  the  CALE  Griineisen  form  violates  Mode  II 
stability.  The  A  =  0  (HULL)  curve  shows  the  DYNA/CALE  curve  for  A  =  0,  which  reduces  into  the 
form  used  by  HULL.  This  form  remains  stable  out  to  the  limit  compression,  p,  =  3,  but  does  not 
follow  the  CALE  curve  well  at  lesser  compressions,  where  experimental  data  were  used  to  fit  the  A 
parameter.  The  P  =  .72  curve,  which  follows  the  functional  form  of  the  currently  proposed  model, 
follows  the  CALE  curve  well  at  lower  compressions,  and  retains  stability  out  to  the  limit  compression. 
Notice  that  the  minimum  value  of  p,  which  is  guaranteed  to  retain  stability,  can  be  obtained  by 
substituting  F0  =  2.09  and  p,  =  1/(1.33  -  1)  =  3  into  the  equation  above.  This  minimum  value  is 
computed  as  0.712,  which  thus  guided  the  selection  of  p  =  .72  for  the  figure. 
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Figure  2.  Comparison  of  GrfingjgenParamcter  Fits  From  CALE  and  HULL  Codes  for  Aluminum 
With  fre  Cumgnqy  Prppoged  Fitting  Model,  in  Liaffl  of  Mode  n  Stability  Criterion. 
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Note  that  some  nen-posilive  values  of  P  arc  pcrmissablc,  but  only  if  T0  <  2/p,  is  true.  For  a 
polynomial  Hugoniot  fit.  where  p,  is  necessarily  infinite,  p  is  necessarily  positive.  One  can  combine 
the  P  constraint  directly  with  the  proposed  T  relation  to  obtain  the  relation  in  terms  of  ro,  p  and  p: 


r  s 


zr0 

2  +  <rfl-2/pJ)  p  ' 


14.  MODE  III  INSTABILITY  USING  THE  M1E-GRUNEISEN  EQUATION  OF  STATE 

Mode  III  instability  arises,  when  the  square  of  the  local  speed  of  sound  is  computed  as  negative — 
a  situation  which  clearly  bears  no  relation  to  any  real  phenomenon.  In  the  Mic-Griinciscn  EOS,  it  can 
arise,  because  of  the  fact  that  the  EOS  is  linearized  about  some  reference  curve,  which  is  the 
Hugoniot,  for  impact  computations.  Since  the  real  world  is  rarely  linear,  this  idealized  linearization 
becomes  less  and  less  accurate  the  further  one  gets  from  the  reference  curve.  When  one  gets  suitably 
distant  from  the  reference  curve,  the  thermodynamic  data  generated  from  a  linearized  EOS  can  be  not 
only  inaccurate  but  also  in  violation  of  basic  thermodynamic  principles. 

When  the  situation  of  imaginary  sound  speed  occurs,  many  codes  (EPIC3,  for  example)  reset  their 
value  to  zero  and  merrily  continue  on  their  computational  way.  Unfortunately,  the  situation  can  be 
much  more  serious  than  just  computing  an  inaccurate  sonic  velocity.  An  imaginary  sonic  velocity 
implies  that  an  increment  of  iscntropic  compression  on  an  element  will  result  in  a  DECREASE  in 
pressure.  Resetting  the  sonic  velocity  to  zero  docs  uuihing  iu  prevent  the  thermodynamic  state  of  the 
element  from  going  berserk.  It  is  these  sorts  of  situations  which  most  often  result  in  catastrophic 
manifestations  of  instability  over  the  span  of  one  or  several  computational  iterations — ballooning 
elements,  collapsing  elements,  and  wild  energy  fluxes  arc  typical. 

The  local  speed  of  sound  is  defined  as  the  rate  at  which  an  infinitesimal  disturbance  propagates 
through  a  medium.  Expressed  by  the  symbol  "C.  it  is  mathematically  expressed  as: 

C*  *  dp/dp  I, . 
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To  derive  its  functional  form  for  the  Mic-Griinciscn  EOS. 


P  ‘  Ph  =  r  P  <E  *  Eh>  > 


take  the  derivative  with  respect  to  an  increment  of  density  (dp)  along  the  isentropc.  This  produces  the 
following  result: 


dp 

3p 


dp  [  dp 


dp  j 


(E  -  £*)  ~  (Fp) 

dp 


Convening  the  derivatives  into  compression  form,  this  relation  may  be  ultimately  simplified,  with 
energy  terms  eliminated,  to  produce  the  result: 


MODE  III  CRITERION: 


Po  C2 


r  + 1 
1  +  p 


r 


r  dVk  Ph  ~Pn 

—  p  *  -  _*  -  °  >  o  . 
2  d  p  1  +  p 


For  the  special  ease  of  p  -  ph  (i.c.,  along  the  Hugoniot),  this  form  expresses  the  relationship  derived 
for  T  in  Section  XI. 


2  \\idpjdy,  -  pp„C2 _ 

p  [\idph/d\k -(pk-p0)f{  1  +  u)]' 


which  was  subsequently  used  to  derive  the  Mode  II  criterion.  Note,  however,  that  Mode  III  only 
requires  the  sound  speed  to  be  positive  at  all  thermodynamic  states,  whereas  Mode  II  requited  the 
slope  of  the  isentropc  (which  is  directly  related  to  the  sound  speed)  to  be  not  only  positive  but 
greater  than  the  Rayleigh  line  slope.  Thus,  the  Mode  II  criterion  is  more  stringent  than  Mode  III,  but 
it  only  applies  along  a  Hugoniot. 

Thus,  if  the  Mode  III  criterion  is  violated  at  some  state,  which  happens  to  lie  on  the  Hugoniot 
fPh-  P).  then  the  Mode  II  criterion  will  have  already  been  violated  at  that  compression  p.  However, 
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the  full  Mode  III  criterion  gives  the  opportunity  to  study  stability  at  all  states  in  the  (p.  v)  or  (p,  p.) 
plane,  not  just  along  the  Hugoniot. 

Operational  simplifications  may  be  possible  when  applying  the  Mode  III  criterion,  depending  on 
the  actual  choices  of  the  and  r  functions.  As  an  example,  if  one  adopts  the  functional  form  for  the 
Oriinciscn  coefficient  proposed  in  Section  XII,  I'  =  ry(l  +  Pg),  then  the  Griinciscn  derivative  terms 
simplify  nicely: 

I  dT  _  I_  n  _  P 

r  dn  r,  (l  +  Pii)  ‘ 


Because  of  the  linearized  nature  of  Mic-GriJneisen,  there  can  be  states  where  Mode  III  stability  is 
violated.  However,  their  occurancc  can  be  greatly  reduced,  given  an  appropriate  choice  of  ph  and  F 
functions.  As  an  example,  consider  Figure  3,  which  depicts  the  (p,  v)  plane  cubic  fit  Hugoniot  for 
aluminum  with  coefficients  .79903,  1.13927,  and  1.39792  Mbar,  respectively,  and  density  of  2.7  g/cm3 
as  given  in  Kohn  (1969).  If  T  is  held  constant  at  its  initial  value  of  2.09,  corresponding  to  p  -  0  in 
the  currently  proposed  model,  the  Mode  III  stability  line,  which  divides  the  (p,  v)  plane  into  stable  and 
unstable  regions,  is  seen  to  produce  instability  in  a  large  part  of  the  (p,  v)  plane,  including 
thermodynamic  states  in  which  the  simulation  might  reasonably  expect  to  exist.  Note  that  an 
isentropc,  when  it  crosses  into  the  unstable  region  of  the  (p,  v)  plane,  changes  its  slope  so  that 
increased  compressions  actually  cause  a  decrease  in  pressure.  It  is  this  sort  of  unstable  behavior  which 
will  cause  a  simulation  to  come  screeching  to  a  hall  in  no  time  at  all. 

On  the  other  hand,  if  P  is  increased  to  a  value  which  guarantees  Mode  II  stability,  as  discussed  in 
the  previous  section,  Mode  III  stability  is  also  greatly  enhanced.  Since  a  cubic  fit  Hugoniot  has  no 
lim-ting  value  of  compression,  the  value  of  P  to  guarantee  Mode  II  stability  must  equal  or  exceed 
TJ2,  which  for  the  ease  in  point  (T0  =  2.09),  is  1.045.  With  this  new  value  of  p,  shown  in  Figure  4, 
the  unstable  region  is  limited  to  a  very  small  portion  of  the  (p,  v)  plane,  which  represents  material 
under  very  large  tensions.  In  fact,  the  tensions  required  to  produce  instability  are  so  large  (391  kbar 
at  v  -  v0,  for  the  case  in  point),  that  aluminum’s  spall  pressure  would  be  exceeded  and  material  failure 
would  occur  before  the  material  was  ever  able  to  approach  this  unstable  thermodynamic  state. 
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Hugoniot 


V/V0 


Depiction  of  Aluminum  =  2.7)  Hugoniot.  Isentropc.  and  Mode  III  Stablilitv  Region 
for  Variable  Grtlneisen  Parameter  (B  =  1.045).  Hugoniot  Is  a  Cubic  Fit.  With  Parameters 
.79903,  1.13927.  and  1.39792  Mbar.  Respectively  '  ’  " 


15.  OTHER  PROBLEMS  WITH  MIE  GRUNEISEN  EOS  IMPLEMENTATIONS 


Movement  of  nodes  along  Lagrangian  coniaci  surfaces  may  produce  occasional  values  of 
compression  in  excess  of  that  anticipated  from  analytical  considerations.  Though  such  artifacts  arc 
purely  computational,  they  can  cause  the  compression  to  go  out  of  (p,  v)  domain  of  known  data  or, 
worse  yet,  out  of  the  physically  plausible  (p,  v)  domain.  When  this  happens,  the  code  may  die 
instantly  or  generate  state  values  so  inaccurate  as  to  cause  computational  "ripples"  which  disturb  the 
rest  of  the  otherwise  valid  computation.  This  problem  can  usually  be  circumvented  by  forcing  a 
smaller  integration  timestep  on  the  computation  or,  alternately,  by  using  an  EOS  formulation  which 
remains  physical  out  to  larger  values  of  compression. 


16.  CONCLUSIONS 


A  review  of  fundamental  shock- transition  theory  was  presented.  The  equation  of  state  was 
introduced  as  an  analytical  vehicle  to  express  material  pressure  as  a  function  of  density  and 
temperature  (or  internal  energy  in  the  ease  of  adiabatic  transition).  The  importance  of  thermodynamic 
principles  was  emphasized  as  a  tool  to  study  the  stability  of  various  EOS  implementations.  The 
Mic-Griincisen  EOS,  with  a  Hugoniot  reference  function,  was  singled  out  for  study  because  of  its 
relative  importance  in  the  computational  modeling  of  shock  transition. 

Starting  from  the  fundamental  laws  of  thermodynamics,  three  criteria  were  developed  to  measure 
the  stability  of  equations  of  state  and  were  applied  specifically  to  the  Mie-Griineisen  EOS,  with 
Hugoniot  reference,  to  investigate  the  stability  characteristics  thereof.  Results  indicate  numerous 
possibilities  of  iastability  under  various  circumstances.  Circumstances  which  could  bring  on  the 
investigated  instabilities  included:  1)  an  improperly  formulated  EOS;  2)  the  choice  of  improper  data 
for  an  otherwise  suitable  EOS;  3)  the  application  of  given  EOS  data  beyond  the  thermodynamic  region 
for  which  the  data  were  originally  intended;  and  4)  the  assumption  of  parameter  constancy  when 
additional  data  were  not  readily  available.  Several  sources  of  published  material  data  and  Hugoniot 
forms  were  investigated  for  stability.  Instability  was  observed  in  these  data,  especially  so  under 
conditions  where  the  applied  compression  was  greater  than  that  for  which  the  data  were  fit 
Unfortunately,  the  existence  in  codes  of  "EOS  material  libraries"  virtually  guarantees  the  indiscriminate 
use  of  EOS  data,  thus  enhancing  the  likelihood  of  EOS  instability. 
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Based  on  this  report’s  findings  and  the  author’s  experience,  it  is  believed  that  a  significant 
percentage  of  hydrocode  computations  which  fail  do  so  as  a  direct  or  indirect  result  of  the  EOS 
calculation.  Typical  code  failures  (for  example:  negative  absolute  energies  or  pressures,  exploding  or 
collapsing  elements,  and  negative  clement  stiffnesses)  arc  due,  probably  in  part,  if  not  completely,  to 
internal  inconsistencies  arising  from  a  misapplied  EOS.  A  basic  understanding  of  thermodynamics  ^nd 
its  relationship  to  shock  transition  will  permit  the  thoughtful  code  developer  and  user  to  verify,  in 
advance  of  computational  solution,  the  validity  of  his  data  and  EOS  forms. 
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Ames  Research  Center 
Moffett  Field,  CA  94035-1000 


10  Dir,  BRL 

ATTN:  SLCBR-DD-T 
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2  Director 
DARPA 

ATTN:  J.  Richardson 

MAJ  R.  Lundberg 
1400  Wilson  Blvd. 

Arlington,  VA  22209-2308 

1  Defense  Nuclear  Agency 
ATTN:  MAJ  James  Lyon 
6801  Telegraph  Rd. 

Alexandria,  VA  22192 

1  Commander 

US  Army  Strategic  Defense  Command 
ATTN:  CSSD-H-LL,  Tim  Cowles 
Huntsville,  AL  35807-3801 

1  Commander 
USA  ARMC 

ATTN:  ATSB-CD,  Dale  Stewart 
Ft.  Knox,  KY  40121 

1  Commander 

US  Army  MICOM 

ATTN:  AMSMI-RD-TE-F,  Matt  H.  Triplett 
Redstone  Arsenal,  AL  35898-5250 

2  Commander 
TACOM  RD&E  Center 

ATTN:  AMCPM-ABMS-SA,  John  Rowe 
AMSTA-RSS,  K.  D.  Bishnoi 
Warren,  Ml  48397-5000 

2  Commander 

US  Army,  ARDEC 

ATTN:  SMCAR-CCH-V,  M.  D.  Nicolich 
SMCAR-FSA-E,  W.  P.  Dunn 
Picatinny  Arsenal,  NJ  07806-5000 

4  Commander 

US  Army  Belvoir  RD&E  Center 
ATTN:  STRBE-NAE,  Bryan  Westlich 
STRBE-JMC,  Terilee  Hanshaw 
STRBE-NAN,  Steven  G.  Bishop 
STRBE-NAN,  Josh  Williams 
Ft.  Belvoir,  VA  22060-5166 
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1  USMC/MCRDAC/PM  Grounds  Wpns.  Br 
ATTN:  Dan  Haywood 

Firepower  Div. 

QuantiCO,  VA  22134 

3  Commander 

Naval  Weapons  Center 

ATTN:  Tucker  T.  Yee  (Code  3263) 

Don  Thompson  (Code  3268) 

W.  J.  McCarter  (Code  6214) 
China  lake,  CA  93555 

2  Commander 

Naval  Weapons  Support  Center 
ATTN:  John  D.  Barber 
Sung  Y.  Kim 
Code  2024 

Crane,  IN  47522-5020 

3  Commander 

Naval  Surface  Warfare  Center 
ATTN:  Charles  R.  Garnett  (Code  G-22) 
Linda  F.  Williams  (Code  G-33) 
Mary  Jane  Sill  (Code  H-11) 
Dahlgren,  VA  22448-5000 

11  Commander 

Naval  Surface  Warfare  Center 
ATTN:  Pao  C.  Huang  (G-402) 

Bryan  A  Baudler  (R-12) 

Robert  H.  Moffett  (R-12) 

Robert  Garrett  (R-12) 

Thomas  L.  Jungling  (R-32) 
Richard  Caminity  (U-43) 

John  P.  Maira 
Paula  Walter 
Lisa  Mensi 
Kenneth  Kiddy 
F.  J.  Zerilli 

10901  New  Hampshire  Ave. 

Silver  Spring,  MD  20903-5000 

1  Director 

Naval  Civil  Engr.  Lab. 

ATTN:  Joel  Young  (Code  L-56) 

Port  Hueneme,  CA  93043 


42 


No.  of 

Copies  Organization 


No.  of 

Copies  Organization 


4  Air  Force  Armament  Laboratory 
ATTN:  AFATUDLJW  (W.  Cook) 
AFATL/DLJW  (M.  Nixon) 
AFATL/MNW  (LT  Donald  Lorey) 
AFATUMNW  (Richard  D.  Guba) 
Eglin  AF8,  FL  32542 

8  Director 

Sandia  National  Laboratories 
ATTN:  Robert  O.  Neliums  (Div.  9122) 
Jim  Hickerson  (Div.  9122) 

Marlin  Kipp  (Div.  1533) 

Allen  Robinson  (Div.  1533) 

Wm.  J.  Andrzejewski  (Div.  2512) 
Don  Marchi  (Div.  2512) 

R  G'aham  (Div.  1551) 

R  Lafarge  (Div.  1551) 

P.  O.  Box  5800 
Albuquerque,  NM  87185 

8  Director 

Los  Alamos  National  l  abcratory 
ATTN  G.  E.  Curt  (MS  K574) 

Tony  Rollett  (MS  K574) 

Mike  Burkett  (MS  K574) 

Robert  Karpp  (MS  P940) 

Rudy  Henninger  (MS  K557,  N-6) 
Roy  Greiner  (MS-G740) 

James  P.  Ritchie  (B214,  T-14) 
John  Bolstad  (MS  G787) 

P.  O.  Box  1663 

Los  Alamos,  NM  87545 

13  Director 

Lawrence  Livermore  National  Laboratory 
ATTN:  Barry  R.  Bowman  (L-122) 

Ward  Dixon  (L-122) 

Raymond  Pierce  (L-122) 

Russell  Rosinsky  (L-122) 

Owen  J.  Alford  (L-122) 

Diana  Stewart  (L-122) 

Tony  Vidlak  (L-122) 

Albert  Holt  (L-290) 

John  E.  Reaugh  (L-290) 

David  Wood  (L-352) 

Robed  M.  Kuklo  (L-874) 

Thomas  McAbee  (MS-35) 
Michael  J.  Murphy 
P.  O.  Box  808 
Livermore,  CA  94550 


1  Battelle  Northwest 

ATTN:  John  B.  Brown,  Jr. 

MSIN  3  K5-22 
P.  O.  Box  999 
Richland.  WA  99352 

1  Advanced  Technology,  Inc. 

ATTN:  John  Adams 
P.  O.  Box  125 
Dahlgren,  VA  22448  0125 

1  Explosive  Technology 

ATTN:  Michael  L.  Knaebel 
P,  O.  Box  KK 
Fairfield,  CA  94533 

1  Rockwell  Missile  Systems  Div. 

ATTN:  Terry  Neuhart 
1800  Satellite  Blvd. 

Duluth,  GA  30136 

1  Rockwell  Intl./Rocketdyne  Div. 

ATTN:  James  Moldenhauer 
6633  Canoga  Ave  rHB  23) 

Canoga  Park,  CA  91303 

2  McDonnell  Douglas  Helicopter 
ATTN:  Loren  R.  Bird 

Lawrence  A.  Mason 
5000  E.  McDowell  Rd.  (MS  543-D216) 
Mesa,  AZ  85205 

1  University  of  Colorado 

Campus  Box  431  (NNT  3-41) 

ATTN:  Timothy  Maclay 
Boulder,  CO  80309 

1  New  Mexico  Inst.  Mining  &  Tech. 
Campus  Station  (TERA  Group) 

ATTN:  David  J.  Chavez 
Socorro,  NM  87801 

2  Schlumberger  Perforating  &  Test 
ATTN:  Manuel  T.  Gonzalez 

Dan  Market 

P.  O  Box  1530/14910  Ariline  Rd. 
Rosharon,  TX  77583-1590 
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2  Aerojet  Ordriance/Exp.  Tech.  C-tr. 

ATTN:  Patrick  Wolf 

Gregg  Padgett 
1100  Bulloch  Blvd. 

Socorro,  NM  87801 

2  Physics  International 
ATTN:  Ron  Funston 

Lamont  Garnett 

2700  Merced  St./P.  O.  Box  5010 
San  Leandro,  CA  94577 

2  Lockheed  Missile  &  Space  Co.,  Inc. 
ATTN:  S.  Kusumi  (0-81-11,  Bldg.  157) 
Jack  Philips  (0-54-50) 

P.  O.  Box  3504 
Sunnyvale,  CA  94088 

i  Lockheed  Missile  &  Space  Co.,  Inc. 

ATTN  Richard  A  Hoffman 
Santa  Cruz  Fac./Empire  Grade  Rd. 

Santa  Cruz,  CA  95060 

1  Boeing  Corporation 

ATTN  Thomas  M.  Murray  (MS-84-84) 

P.  O  Box  3999 
Seattle.  WA  98124 

2  Mason  &  Hanger  -  Silas  Mason  Co. 
ATTN:  Thomas  J.  Rowan 

Christopher  Vogt 
Iowa  Army  Ammunition  Plant 
Middletown,  IA  52638-9701 

1  Nuclear  Metals  Inc. 

ATTN:  Jeff  Schreiber 
2229  Main  St. 

Concord,  MA  01742 

1  Lockheed  Engineering  &  Space  Sciences 
ATTN:  Ed  Cykowski,  MS  B-22 

2400  NASA  Road  1 
Houston,  TX 

2  Dyna  East  Corporation 
ATTN:  P  C.  Chou 

R.  Ciccarelli 
3201  Arch  St. 

Philadelphia,  PA  19104 


2  Southwest  Research  Institute 
ATTN:  C.  Anderson 
A.  Wenzel 
6220  Culebra  Road 
P  O.  Drawer  28510 
San  Antonio,  TX  78294 

2  Battelle  -  Columbus  Laboratories 
ATTN:  R.  Jameson 

S.  Golaski 
505  King  Avenue 
Columbus,  OH  43201 

3  AUiant  Techsystems,  Inc. 

ATTN:  Gordon  R.  Johnson 

Tim  Holmquist 
Kuo  Chang 
MN  48-2700 
7225  Northland  Dr. 

Brooklyn  Park,  MN  55428 

1  S-Cubed 

ATTN:  R  Sedgwick 

P.O.  Box  1620 

La  Jolla,  CA  92038-1620 

2  California  Research  and  Technology, 

Inc. 

ATTN:  Roland  Franzen 
Dennis  Orphal 
5117  Johnson  Dr 
Pleasanton,  CA  94566 

2  Orlando  Technology,  Inc. 

ATTN:  Dan  Matuska 

J.  Osborn 
P.  O.  Box  855 
Shalimar,  FL  32579 

3  Kaman  Sciences  Corporation 
ATTN:  D.  Barnette 

D.  Elder 
P.  Russell 
P.  O.  Box  7463 
Colorado  Spring,  CO  80933 
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2  Defense  Research  Establishment  Suffield 
ATTN:  Chris  Weickert 

David  MacKay 

Ralston,  Alberta,  TOJ  2NO  Ralston 
CANADA 

1  Defense  Research  Establishment  Valcartier 
ATTN:  Norbert  Gass 
P.  0.  Box  8800 
Courcelette,  PQ,  GOA  1RO 
CANADA 

1  Canadian  Arsenals,  LTD 
ATTN:  Pierre  Pelletier 
5  Montee  des  Arsenaux 
Villie  de  Gardeur,  PQ,  J5Z2 
CANADA 

1  Ernst  Mach  Insiitute 
ATTN:  A.  J.  Stilp 
Fckerstrasse  4 
D-7000  Freiburg  i.  Br. 

GERMANY 

3  IABG 

ATTN:  H  J.  Raatschen 
W.  Schittke 
F.  Scharppf 
Einsteinstrasse  20 
D-8012  Otlobrun  B.  Munchen 
GERMANY 

1  Royal  Armament  R&D  Establishment 
ATTN:  Ian  Cullis 
Fort  Halstead 

Severioaks,  Kent  TNi4  7BJ 
ENGLAND 

1  Centre  d'Etudes  de  Gramat 
ATTN:  SOLVE  Gerard 
46500  Gramat 
FRANCE 

1  Centre  d'Etudes  de  Vaujours 
ATTN:  PLOTARD  Jean-Paul 
Boite  Postale  No.  7 
77181  Country 
FRANCE 
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1  PRB  S.A. 

ATTN:  M.  Vansnick 

Avenue  de  Tervu6ren  168,  Bte.  7 

Brussels,  B-1150 

BELGIUM 

1  AB  Bofors/Ammunition  Division 
ATTN:  Jan  Hasslid 
BOX  900 
S-691  80  Bofors 
SWEDEN 
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